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A model is proposed where two x' 2 ' nonlinear waveguides are contained in a cavity suited for 
second-harmonic generation. The evanescent wave coupling between the waveguides is considered 
as weak, and the interplay between this coupling and the nonlinear interaction within the waveg- 
uides gives rise to quantum violations of the classical limit. These violations are particularly strong 
when two instabilities are competing, where twin-beam behavior is found as almost complete noise 
suppression in the difference of the fundamental intensities. Moreover, close to bistable transitions 
perfect twin-beam correlations are seen in the sum of the fundamental intensities, and also the self- 
pulsing instability as well as the transition from symmetric to asymmetric states display nonclassical 
twin-beam correlations of both fundamental and second-harmonic intensities. The results are based 
on the full quantum Langevin equations derived from the Hamiltonian and including cavity damp- 
ing effects. The intensity correlations of the output fields are calculated semi-analytically using a 
linearized version of the Langevin equations derived through the positive-P representation. Con- 
firmation of the analytical results are obtained by numerical simulations of the nonlinear Langevin 
equations derived using the truncated Wigner representation. 
PACS: 42.50.Dv, 42.50.Lc, 42.65. Sf, 42.65.Wi 



I. INTRODUCTION 

The nonlinear materials have been the subject of 
various investigations in recent years. Using a cavity 
setup the weak nonlinearities can be resonantly ampli- 
fied, and complex spatiotemporal behavior has been pre- 
dicted from a classical point of view [Q. Moreover, due to 
the quantum fluctuations of light many interesting non- 
classical effects have been reported, such as squeezed light 
and sub-Poissonian light 0, both theoretically 
and experimentally |j,|7j . The interplay between the clas- 
sical spatial instabilities and the quantum fluctuations in 
the system has been investigated intensively lately |l| , a 
study devoted to characterizing the mode interaction on 
the quantum level. 

We consider the case of second-harmonic generation 
(SHG), where the photons of the pump field (funda- 
mental photons) are up-converted in pairs to second- 
harmonic photons of the double frequency. The model 
we propose in this paper consists of two quadratically 
nonlinear waveguides placed in a cavity that resonates 
both the fundamental and second-harmonic, and we take 
linear coupling between the waveguides into account. In 
a sense this is the simplest mode coupling model ob- 
tainable. The question is how the coupling between the 
waveguides affects the cavity dynamics, and in particular 
we shall focus on the nonclassical behavior of the system. 

The name proposed for this model, the quantum opti- 
cal dimer, originates from the numerous investigations 
made about discrete site coupling in various systems, 
such as condensed matter physics and biology, see Ref. 
H for a general treatment of discrete systems. Thus, the 
name dimer implies that coupling between two discrete 



sites are being taken into account. For a single waveguide 
(or site) we shall use the name monomer, a case corre- 
sponding here to a bulk nonlinear medium and neglecting 
diffraction. 

It has been shown that in the SHG quantum optical 
monomer excellent squeezing of the output fields is pos- 
sible close to a self-pulsing instability H, and nonclassi- 
cal effects in SHG have been verified experimentally [Q, 
and even shown to persist above the threshold of the in- 
stability Also in the presence of diffraction strong 
correlations exist between different spatial modes in the 
presence of a spatial instability jlCj , including strong cor- 
relations between the fundamental field and the gener- 
ated second-harmonic field as well as spatial multimode 
nonclassical light. 

As we shall show the quantum optical dimer also dis- 
plays strong nonclassical correlations in the intensity cor- 
relations, and that the linear coupling across the waveg- 
uides plays a decisive role. The model has three types of 
transitions, bistability, self-pulsing and a symmetric to 
asymmetric transition. It is remarkable that particularly 
strong nonclassical correlations are observed when two 
of these instabilities compete. Specifically, when taken 
close to a self-pulsing or bistable regime the symmetric 
to asymmetric transition has nearly perfect twin-beam 
behavior, so the difference of the fundamental intensities 
displays almost no fluctuations. The twin-beam correla- 
tions were first shown in the optical parametric oscillator 
(OPO) pl] , P^ i where the signal and idler photons of the 
twin-fields are generated simultaneously from the pump 
field, and the intensity difference shows correlations be- 
low the classical limit. However, the twin-beam effect 
observed in the dimer originates from photons created 
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in different waveguides with only the coupling to link 
them. Thus, the photons are strictly speaking not twins, 
but merely "brothers" . The twin-beam effect is also ob- 
served near bistable turning points where complete noise 
suppression is observed in the sum of the fundamental 
intensities, the strongest violations occurring in the limit 
where the fundamental input coupling loss rate is much 
larger than the second-harmonic one. That bistability 
gives rise to highly nonclassical effects turns out to also 
hold for the SHG monomer, and has to our knowledge not 
been observed before; usually the self-pulsing transition 
has been used to observe violations of the classical limit, 
which we also observe in the dimer. The bistable transi- 
tion has previously been observed to produce nonclassical 
states in other systems, such as dispersive and absorptive 
optical bistability and Raman lasers [Q . 

A closely related optical model is spatially coupled 
lasers ]f5| , where a single laser medium is pumped by 
two beams spatially separated. Waveguiding is achieved 
by thermal lensing, in which the temperature dependent 
refractive index of the laser medium creates a guiding 
effect, and the coupling strength is controlled by the dis- 
tance between the pump profiles. The quantum noise 
induced correlations in these systems have not yet been 
reported to beat the standard quantum limit when Kerr 
type nonlinearities are considered jL&jl7|, except when 
the coupling arises solely due to initially correlated noise 
terms of the pumps |jq| . 

The cavityless setup of coupled waveguides has 
previously been investigated, both from a classical and 
a quantum mechanical point of view. In the classical 
model of waveguide arrays, the focus of attention has 
been on soliton behavior originating from the coupling 
fl9fl , whereas the cavityless dimer was shown to pro- 
duce chaotic states away from the integrable limit (where 
second-harmonic coupling is neglected) ||(|. The quan- 
tum behavior of the cavityless dimer has been investi- 
gated by the group of Perina et al. (for a review see Ref. 
]2l|) giving the name "nonlinear coupler" to the model. 
They have investigated both co- and counter-propagating 
input fields in parametric oscillation and, e.g., the trans- 
fer of quantum states from one waveguide to the other. 

The model presented here is also closely related to the 
dynamics of coupled atomic and molecular Bose-Einstein 
condensates (BECs) |2^,|2j|, where the photoassociation 
of an atomic condensate may produce a molecular con- 
densate with an atom-molecule interaction that is remi- 
niscent of the interaction between the fundamental and 
second-harmonic photons in SHG J24J . The opposite pro- 
cess where the photo dissociation of a molecular BEC cre- 
ates an atomic BEC has been shown to produce squeezed 
states ]2q| , a model which has the quantum optical equiv- 
alent in the OPO. If an analogy should be drawn between 
the quantum optical dimer presented here and BEC it 
would consist of placing two such coupled molecular- 
atomic BECs in separate quantum wells. Thus, evanes- 
cent tunneling of the wave functions between the wells 
would introduce the dimer coupling, similar to what is 



done in Ref. [g6|] for a normal BEC. 

We should finally stress that the cavity setup discussed 
in the present work gives rise to two major differences 
to the work in cavityless waveguides as well as for the 
BEC. First of all the cavity introduces losses in the model 
through the input mirror, and secondly external pump 
fields appear in the equations acting as forcing terms. 

The paper is structured as follows. In Sec. ||the model 
is introduced, and the stochastic Langevin equations are 
derived from the full boson operator Hamiltonian. Also, 
we discuss the allowed values of the coupling constants. 
In Sec. Ill the linear stability of the Langevin equa- 
tions are investigated, and the bifurcation scenario of the 
model is discussed. In Sec. [V we discuss the framework 



for the two-time photon number correlations of the out- 
put fields, and the semi-analytical spectral variances are 
derived in the linearized limit. Section^ is devoted to the 
results of the analytical calculations as well as the numer- 
ical simulations. A summary is made in Sec. VI where 
we also discuss the results obtained. Appendix~|A shows 
details about the derivation of the quasi-probability dis- 
tribution equations used to connect the master equation 
for the quantum Hamiltonian with the classical looking 
stochastic Langevin equations. The numerical methods 
are discussed in App. H. 




FIG. 1. The setup. Two nonlinear waveguides A and B 
inside a cavity pumped by a classical field. 



II. THE MODEL 



We consider the setup shown in Fig. [|. Two \^ non- 
linear waveguides are contained in a cavity with a high 
reflection input mirror M\ and a fully reflecting mirror 
Mi at the other end. The cavity is pumped at the fre- 
quency u>i and through the nonlinear interaction in the 
waveguides photons of the frequency u>2 = 2wi are gener- 
ated; this is the process of second-harmonic generation. 
The cavity supports a discrete number of longitudinal 
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modes, and we will consider the case where only two of 
these modes are relevant, namely the mode wi jCav closest 
to the fundamental-harmonic (FH) frequency and W2,cav 
closest to the second-harmonic (SH) frequency. Using 
the mean field approximation the z-dircction, in which 
the pump beam propagates, is averaged out. This ap- 
proach is justified as long as the losses and detunings are 
small, and we furthermore assume perfect phase match- 
ing. The waveguiding implies that diffraction in the 
transverse plane may be neglected. Let A\(t) and Bi(t) 
(A 2 (t) and -82(4)) denote the FH (SH) intracavity boson 
operators of waveguide A and B, respectively. They are 
normalized so they obey the following equal time com- 
mutation relations 

[d i (t),6}(t)]=s ij , i,j = 1,2, 6 = a,b, (i) 

while [Aj(t),Bj(t)} = 0. The syste m is modelled through 
the Hamiltonian 



H 



Hab, 



(2) 



where the system Hamiltonians in the frame rotating 
with the pump frequency are given by 

h s J s = -Mi 61 61 - m 2 6\6 2 + l -^(6\ 2 6 2 - 6161) 



-ih{£ p ,oO\-£; 0{), = A,B. 



(3) 



The detunings from the cavity resonances are given by 
Sj = ujj — Wj lCa v, K is proportional to the \^ nonlin- 
earity and £ p ,o are the external pump fields at the FH 
frequency of the individual waveguides |^7j , here treated 
as classical fields. The coupling between the waveguides 
is modelled as overlapping tails of evanescent waves so it 
may be assumed weak, implying we can describe it as a 
linear process 



Hab = hJ^AxBl + + hJ 2 {A 2 B,\ 



B2M) 



Ji and J 2 are the cross-waveguide coupling parameters 
of the FH and SH, respectively. The time evolution 
of the reduced system density matrix operator p in the 
Schrodinger picture is then given by the master equation 

MM 



dp 
m 



:[H,p] + {Lx.a + L 2 ,a + Li.b + L 2 . B )p. 



(5) 



The continuum of modes outside the cavity is modelled 
as a heat bath in thermal equilibrium, and the coupling 
to these modes has been included through the Liouvillian 
terms 

+7^f([6 i A6j] + [6],y36 i ])- (6) 



These terms describe the losses of the fields through pho- 
tons escaping the cavity, and simultaneously they model 
fluctuations entering the cavity through the input mir- 
ror, a consequence of the dissipation-fluctuation theorem 
p0| . The loss rates of the input coupling mirror are given 



by 7j, whereas the terms n* h 



^ e huij/k E 



1) 



the mean number of thermal quanta in the external bath 
modes at ujj. We shall here neglect thermal fluctuations 
by setting the bath temperature T — yielding n* h = 0. 
First of all this is a good approximation for optical sys- 
tems since here hui > fcgT, and secondly we may hereby 
focus on behavior solely due to the inherent quantum 
fluctuations of light. 

The master equation (|J) is difficult to solve as it is, 
therefore we apply the now standard technique of ex- 
panding the density matrix in a basis of coherent states 
weighted by a quasi-probability distribution (QPD). The 
details of this quantum-to-classical description are given 
in App. |A|, and the result is a partial differential equa- 
tion of the QPD. This QPD equation depends on the 
choice of ordering of the c orr esponding quantum mechan- 
ical averages. Equation (AS) is the QPD equation using 
the positive-P distribution giving normally ordered av- 
erages, which we will use for the linearized analysis. For 
the numerical implementation the Wigner distribution is 
used to obtain Eq. (|A^) , in which symmetric averages are 
calculated. 

If the QPD equations (A8) and (|A9) are on the 
Fokker-Planck form ( A10 ), an equivalent set of stochastic 
Langevin equations (All) can be found to by using Ito 
rules of stochastic integration |nj . For the Wigner QPD 
equation ( jASj ) this is not the case because of the third or- 
der terms, however these terms, which have been shown 
to model quantum jump processes |32| , are generally ne- 
glected and the resulting Fokker-Planck equation turns 
out to be a good approximation to the original problem. 
Using this approximation the normalized Langevin equa- 
tions for the Wigner QPD equation are 



( 4 ) A 1 = (-l + iA 1 )A 1 +AlA % -iJ l Bi+V2Ai a ,i(t) (7a) 
A 2 = (-7 + iA 2 )A 2 - \^A\ - iJ 2 B 2 + y^A ina (t) (7b) 



Bi = (-1 + iAi)Bi + B{B 2 - U X A X + s/2B iliA (t) (7c) 
B 2 = (-7 + iA 2 )B 2 - X -B\ - iJ 2 A 2 + y^B in , 2 {t), (7d) 

where the dot denotes derivation with respect to time. 
The fields {Aj, A*} and {Bj, B*} are normalized equiva- 
lent c- numbers to the operators and {Bj,Bj}. 
The input fields are describing the pump field entering 
the cavity through the input mirror as well as the noise 
coupled in here according to the Liouvillian terms (|^) 

F in<1 (t) = JL+£ Fl (t), F in>2 (t) = £f 2 (t) (8a) 

mmFA t ')) = s i /-^^, (8b) 
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with F — A, B. All oth er c orrelations are zero. The 
positive- P QPD equation (|A8[ ) is on Fokker-Planck form 
so no approximations are needed. The equivalent set of 
Langevin equations is given by (Q) with A* — > Aj and 

B* —¥ Bj, as well as the equations for the fields 



A\ = (-l-iAi)Al 

+i.hB\ + V2A\ al {t) 

4 = (- 7 -iA a )4-i(4) a 

+ij 2 i?| + v^4 ni2 (i) 

Si = (-1 - iAi)_B| + 

+iJiA\ + V2Bl nl (t) 

^-(-J-iA 2 )Bl-l(Blr 

+U 2 Al + ^BtJt). 



S 2 = 



(9a) 

(9b) 
(9c) 

(9d) 



The fields {A V A]} and {B^B]} are normalized equiva- 
lent c-numbers to the operators {Aj,Aj} and {Bj,Bj}. 
The input fields for the positive-P Langevin equations 
are 

Fin,i(t) = -j= + fa(t), <i(*) = ^+4W ( 10a ) 
Pin, 2 (i)=0, < 2 (t)=0 (10b) 

(^(^(0) = ^*"° (10c) 

(4(^)4(0) = ^?''^ ( 10d ) 

with F — A,B, and again all other correlations are zero. 
The doubling of phase space associated with the positive- 
P representation (see App. |X] for details) implies that 
£] is uncorrelated to and also that Ai and A\ are 
independent complex numbers and only on average is 
A]=A*. 

The Langevin equations have been normalized by in- 
troducing the dimensionless variables 

i = 7l t, 7 =^, A, = ^ (11a) 
7i 7i 



-4, 



7i 



Ai n> j'(t) 



3/2 

7i 



Q!in,j(*), -Binj(i) 



Bi = -^ (Hb) 

7i 

-^A n ,iW (He) 

7i 



-4 



(lid) 



and the tildes have been dropped. The fields aj and (3j 
are the unsealed c-numbers, cf. App. |A|. We have fur- 
thermore introduced the dimensionless quantity 



This parameter sets the level of the quantum noise, cf. 
Eqs. (§) and ©, and in the OPO it represents the sat- 
uration photon number to trigger the parametric oscilla- 
tion. 

For simplicity we have assumed real and equal pump 
rates in both waveguides £ Pt A = £p.B = £p- The con- 
sequence is that the same input mirrors as well as in- 
tracavity paths are used for both waveguides, implying 
identical detunings Q as well as losses for the FH fields 
and the SH fields, respectively. 
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FIG. 2. Parallel planar waveguide setup shown in the 
:r-direction. The width of the waveguides is Z w and the dis- 
tance between them is d. The transverse distributions of the 
lowest order modes for a realistic setup are shown calculated 
using a step profile of the refractive index. 

The coupling strengths between the waveguides are 
controlled by J\ and J 2 , and it is relevant to consider 
what values these may take. Fig. || shows an instructive 
example, where we consider symmetric step-index par- 
allel planar waveguides with a core (cladding) refractive 
index n co (n c \). The weakly guiding limit is assumed 
where n co ~ n c \ . Taking the FH field of waveguide A as 
example, the coupling from waveguide B can be found 
by considering waveguide A in isolation and taking the 
presence of waveguide B as a weak perturbation. This 
approach assumes that the transverse profile u(x) and 
propagation constant [3 of the modes in the waveguides 
are left unchanged, and only the amplitude is modified by 
the perturbation. The coupling constants of the propa- 
gation equations of the waveguides are then found as p4|] 



Ja x Bi — 



t 2 k 2 

l cl ft l 



-d/2 



dx u Al {x)u Bl {x), 
PAi J-d/2-U- (13) 



= -77?. 



(12) 



where k\ — 2tt/\i is the vacuum wavenumber of the 
FH, and the mode profiles are assumed normalized so 
dxu 2 (x) = 1. Thus, Eq. ( fl3| ) has the dimension 

m _1 . Fig. ^ shows the lowest order modes of the isolated 
waveguides as calculated for a realistic setup. If only 
coupling between modes of same order is considered we 
have Ja 1 b j _ = Jb^Ai = Jf r ° p - Applying the mean field 
approach Q the coupling parameters of Eqs. (]?]) is then 
given by Jj — Jj mp L ca _ v /Ti, where L cav is the length of 
the c avit y, and r is the cavity round trip time. From 
Eq. (lid) the normalized coupling parameter is found 
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through 7x = T\/(2t) where T\ is the FH intensity trans- 
mission efficiency of the input mirror, so we obtain |36| 



J, = J: 



prop 2i C av 



(14) 



As a result of these considerations we see that the SH 
coupling parameter generally will be lower than the FH 
one. This is clear from the calculated modes in Fig. 
where the SH modes (dashed) decay faster than the FH 
modes (solid). However, it is impossible to generally say 
how much weaker and when the distance between the 
waveguides is decreased the coupling parameters become 
closer to each other. Finally, the actual values of Jj are 
highly sensitive to the specific setup. Not only in terms of 
waveguide parameters (e.g., distance between guides, the 
modes in the guides), but also on independent parame- 
ters (cavity length, input transmission efficiency). Using 
parameters from realistic setups (similar to the cavity 
setup discussed in Ref. p7|) we easily obtained normal- 
ized coupling parameters of O(10 2 ), while still preserving 
the assumptions of weak coupling as well as the mean- 
field limit. Finally, when coupling between the lowest 
order modes is considered we have Jj > 0. 



III. LINEAR STABILITY 



In this section the Langevin equations derived previ- 
ously are linearized and the linear stability is investi- 
gated to obtain a bifurcation scenario in the classical 
limit where noise is absent (n s — > oo). In this limit 
the Langevin equations from the different representations 
give the same result, a natural consequence from the fact 
that in the classical limit the operators commute. Ad- 
ditionally, in Sec. |y] we are going to use the linearized 
equations with noise to derive analytical results for the 
noise induced correlations. For this purpose it is more 
convenient to use normally ordered intracavity averages, 
as will be explained later, and this section will therefore 
only concern the positive-P Langevin equations. The 
results of this section reveal both symmetric and asym- 
metric modes in the two waveguides, as well as bistable 
behavior and Hopf unstable solutions. 

The linearization is particularly simple in the symmet- 
ric case. Here the steady states in the waveguides are 
identical, so the FH steady states in waveguide A and B 
are equal and equivalently for the SH steady states. The 
symmetric steady states of the waveguides can be found 
from the monomer equations, i.e. using the results of 



Ref. [B8| and apply the substitution A, 



J, 



In the symmetric case the steady states are denoted 



E z = I 



+ h{d\ + l) 



(15a) 

(15b) 
(15c) 
(15d) 

We may linearize the positive-P Langevin equa- 
tions ([?]), (0) and ( |i~0| ) around the symmetric steady 
states Aj = AAj + A and B = ABj + A [|§ to get 
the matrix equation 



I 2 =P 1 [4(dl + 1 2 )]-' 

01 = -Arg (1 - idx + Ji/[2( 7 - id 2 )]) 

02 = -Arg(- 7 + id 2 ) + 20i. 



Aw = AA\ 



B 



n(t), 



(16) 



where Aw is a vector of fluctuations 



Aw 



A Ai , A A\ , A^ 2 , AAl APi , AB\ , AP 2 , ABl 



and n(i) is a vector of Gaussian white noise terms corre- 
lated as 

(n j (t)n k (t'))=S(t-t'). (17) 
The matrix A is block ordered into four 4x4 matrices 



A = 



A A 

A. T A n 



with the diagonal cross-coupling matrix 

A x = Diag [ —iJ\ %J\ —iJ 2 iJ 2 ] , 
and the monomer matrix is 





"-1 + idi 


A 2 


At 


■ 


A 


At 


— 1 — id\ 





A x 


^-m — 


-A, 





-7 + id 2 










-At 





-7 - id 2 _ 



Aj = Bj = x T,< 



giving 



(18) 



(19) 



(20) 



The diffusion matrix is also diagonal 

D = Diag [A 2 Al A 2 A* 2 0] , (21) 

and D = B T B. 

The classical stability of the system is found by solv- 
ing the eigenvalue problem Av = Av, which is done 
in Mathematica. The analysis is characterized by 
two cases, either w hen o n e ph ysical solution exists to 
the closed problem (15a)-( [15b| ), or when the system is 
bistable and three physical solutions exist (in this case 
each solution must be analyzed individually). The sta- 
bility of the steady states may now either change with the 
critical eigenvalue Aj iC having Im(Aj. c ) = at the critical 
pump value P S ym, which means that the symmetric state 
of the system is no longer stable. When this happens a 
new state with Aj ^ Bj is stable instead, and the actual 
values of these new steady states are not easily calculated. 
We will not address the stability of the system beyond the 
asymmetric transition any further in this paper, however 
the transition to the asymmetric state will be used to 
look for nonclassical correlations. The other possibility 
is that the system changes stability with lm(Xj C ) ^ at 
the critical pump value i?sp , which corresponds to a Hopf 
instability leading to self-pulsing temporal oscillations. 
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The system is well characterized by the relative loss 
rate 7, which in the SHG monomer was shown to de- 
termine the degree of squeezing as well as which field 
the best squeezing was observed ||. Following the sim- 
ple layout of the cavity shown in Fig. [I] implies that 
Ai = A2 = A JIU. The bifurcation scenario in the 
{Ji,A} space for J2 = 1.0 is shown in Fig. ||, which 
displays a rich variety of behavior. For A < 1 bistable 
behavior is observed, and the upper branch may be both 
Hopf unstable as well as asymmetrically unstable as in- 
dicated. For A > 1 a large Hopf region is seen, while for 
Ji large, asymmetric states are observed. 

Setting 7 = 1 a similar scenario as for 7 = 0.1 is ob- 
served: On resonance self-pulsing symmetric states dom- 
inates, while bistable solutions may be seen for A < — 1 
and asymmetric states appear when A > 0. For 7 = 10 
the self-pulsing instability dominates and asymmetric so- 
lutions only appear when detuning is introduced and si- 
multaneous large values of J\ and J 2 are chosen. Bistable 
solutions are not seen here except for very large coupling 
strengths, a consequence of the criteria for bistability 
(Eq. (6) in Ref. M) 



\d 2 \{\di\-V3) 



> 7, did 2 > 0. 



(22) 



V3\dx\ + 1 

All these results indicate that the most diverse bifurca- 
tion scenario is when 7 < 1. 



spectrum of fluctuations in the stationary state, provided 
that the fluctuations are small. These intracavity fluctu- 
ation correlations can be directly related to the output 
correlations by using the input-output theory of Gardiner 
and Collett [^(J . We will only present results in the case 
where the symmetric steady states are stable. 

The input fields 0- m j(i) coupled into the cavity 
through the input mirror are posing an instantaneous 
boundary condition for the output fields 



O, 



out J' 



(*) 



OinjC*), = A,B. 



(23) 



It should be stressed that in this equation 7^ is only the 
loss rate of the input mirror, and does not include addi- 
tional absorption losses of the cavity that might other- 
wise have been included in the Langevin equations. Also 
note that the input operator is in the FH taken as a both 
the classical pump as well as the fluctuations around this 
classical level originating from the heat bath interaction 
[so really an operator equivalent of the Langevin input 
fields from Eqs. (§) and (0)]. The fields outside the 
cavity obey the standard free field commutator relations 



[OontAt),dt ut dt')] = s jk s(t-t'), 



(24) 



and 




10 
J, 



15 



20 



(25) 



while all other commutators are zero. We want to express 
correlations of the out fields entirely on correlations of the 
intra cavity fields, hence we want to get rid of terms in- 
volving Oi n .j(t). Using arguments of causality it may be 
shown that this can only be done if time and normally 
ordered correlations are considered |^(|, as e.g. 

(^LtjW. ^WO) = 2 7j -(it(t),i j (t')) (26a) 
(A outJ {t),A outJ (t')) = 27 i (i j (max[t,t']), 

i^min^t'])), (26b) 



FIG. 3. Bifurcation diagram showing the stability for 
Ai = A2 = A, 7 = 0.1 and J2 = 1.0. In the bistable area the 
stability of the upper branch is indicated. 



IV. PHOTON NUMBER SPECTRA 

The linearized Langevin equation for the positive-P 
representation can be used to analytically calculate the 



which precisely implies time and normal order of the cor- 
relations. Since this is exactly what the P-representation 
computes the intracavity operator averages on the right 
hand side of Eq. ( pt^ ) may be directly replaced by c- 
number averages from the P-representation. 

The intensities of the output beams may be found from 
the photon number operator N^ ut j — Ol ut j0 out j. In a 
photon counting experiment two-time correlations of the 
intensities may be calculated as 



CT> Bk (r) = (N* ttj (t) ± N» tik (t),N* ttj (t + t) ± N^ k (t + r)) 

= ii^Lj) + (N^ t!k ))S(r) + (: N* tJ (t) ± N^ k (t), N* tJ (t + r) ± N* t , k (t + r) :) 



2(7i<#/> +lk(N k B ))5(r) +4(: 6N ( A %(t),5N 



(27) 



G 



where the notation (: :) indicates a normally ordered 
average. In the second line we have used the commutator 
relations (53) to rewrite to normal order, while the last 
line follows from Eq. (|26|). We have also introduced 



V 



(±) 



(*) = -YjN?(t) ± l k N?(t), Nf = ()](),. 



o - nt, 



(28) 



and calculated the variance as (S,T) = (ST) - (S)(T). 

It is more convenient to investigate these two-time 
correlations in the Fourier frequency domain using the 
Wiener-Khintchine theorem p0| 

/oo 
-OO 

= 2( lj (Nf)+ lk (N k B }) 



4 / dre 1 ' 



5N%l h (t),5N%l k (t + T)-.) (29a) 



(29b) 



Here we have introduced the spectrum normalized to 
shot-noise vj^g k (to), and the shot-noise level given by 

C S N=2( lj {Nf)+ lk {N*}), (30) 

is with this normalization unity. The shot-noise level 
is the limit between classical and quantum behavior, 
hence if Aj and B k are coherent states the variance will 
be vj^g k (uj) = 1. A complete violation of the shot- 
noise level Va~b (lo) = implies that no fluctuations 
are associated with the measurement of the intensities 
N^ t j ± N^ nt k . The correlations between the fields of 
the same waveguide are 

= 2 lj {N 3 A )6(r)+i^(:N Aj (t),N Aj (t + r) :), (31) 
which means that the monomer spectra are 



V Aj {u)=2 lj {Nf) 



/ dre WT {: N Aj (t),N Aj (t + r) :), (32) 



so the shot-noise level is here Csn = 2jj(Nf). 

Until now everything has been kept in operator form. 
The next step is connecting the operator averages with 
c-number classical averages, which will here be done with 
the semi-analytical calculations in mind. Thus, we apply 
the positive- P representation averages and note that we 
shall only consider symmetric states making Aj = Bj, 
and that the spectra eventually calculated are linearized. 

Expressing the spectra ( |29a| ) and ([32]) in dimensionless 
c-numbers from the P-representation we readily have 



(uj) = 2n s 1 {j j I j +%I k ) 



dre WT {5I 



kk^'^kk 



(< + r)) P (33a) 



V Aj (w) = 2^1 



dre^(lf(t),lf(t + r)) 



(33b) 



with the subscript P referring to the averages being cal- 
culated in the P-representation. We have here used 
7j = 7j/7i and the c-number equivalent of Eq. ( p8| ) 



si A % k (t) = vif(t)±*f k i*(t), if in-,. 



(34) 



while the shot-noise level (|30|) is 

Csn =2n- 1 (*/ j I j +%h). (35) 
The dimensionless spectra are found by the searings 

(36) 

and tildes have been dropped. 

The linearized equations ( jig) m ay be solved directly 
in frequency space (see Ref. }4l|] for details). So let 
us define the spectral matrix of fluctuations in the P- 
representation in the steady state limit (where we may 
choose the time t arbitrarily and henceforth take t = 0) 



/oo 
dre Mr (Aw(0)Aw(T) T )j 
-OO 



(37) 



with the superscript n indicating that the averages are 
equivalent to normally ordered quantum mechanical av- 
erages. This may be calculated using 

S n (oj) = (-iwl - A) -1 D(jwI - A T )~\ (38) 

where I is the identity diagonal matrix. 

In order to calculate intensity correlations we evaluate 
terms like (1^(0), I k {r))p, and expressing this in terms 
of the fluctuations around the symmetric steady state, 
second, third and fourth order correlations in Aw are 
obtained. Due to the strength of the steady state values 
higher order correlation terms may be neglected, so we 
get to leading order 

(//(0),/ fc s (r))p ~ AjA* k {AA]{Q),AB k {r))p 
+A*A k (AAj-(O), AB ] k (r)) P + A* A* k (A A, (0) , AB k (r)) P 
+A J A k (AA](0),ABl(r)) P . 

U sing this result the normalized dimensionless spec- 
tra (33a) are 
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V£i» = 1 + _-_{y£( w ) + 7 2 V%» ± 2 7 Re[^^l 2 (5 1 n 4 ^) + ^WJ + ^fS&M + S&H)] } (39a) 

- 1 + Tr-^-yrjvy, (") + 7 2 ^, M ± 2 7 Re[AlA 2 (S? 8 (u J ) + S&(-w)) + A\A* 2 {S? 7 {u) + S^(u))] } (39b) 

11 + "fl2 1 J 

= 1 + ^ ; 7 - ^ ± 2Re[5r 6 (o;) + S5,(- W ) + e^^M + (39c) 



± 2 7 Re[^ 8 (u,) + S&(- w ) + e^S^) + S%(w))] . 



(39d) 



Here we have used the normally ordered (indicated 
with a superscript n) single mode spectrum defined as 



V?» = I dre^(lf(0),lf(r)) P , (40) 

V2M 



V2M 



/i(sT 2 (o;) + ^ a (-w) 
+2Re[S*i\(u;)e-^ 1 

7 2 (5 3 " 4 H + 5^ 4 (-c) 
+2Rc[^ 3 (o;)e- 42 ^ 



(41a) 



(41b) 



and Vg . (w) = V^ (lo). With these quantities the 



monomer spectra (33b) normalized to shot-noise are 
readily calculated 



V Aj {u>) = l + ^V2.(u;). 

1 i 



(42) 



The calculations of the spectra use the general sym- 
metry properties of the spectral matrix S n (u>), so e.g. 
S?i (w) = [S 2 " 2 (w)]* and (w) = ^ (-w). 



V. RESULTS 

In this section we present intensity correlation spec- 
tra both from the semi-analytical derivation, as well as 
results from the numerical simulations (the numerical 
method is discussed in App. |b|). The chosen examples 
are only illustrative for the overall behavior, and the re- 
sults hold for large parameter areas. This is especially 
important to stress for the coupling parameters, since 
they are not so easy to control experimentally as the de- 
tunings and loss rates. 

In order to understand the effect of the coupling be- 
tween the waveguides, a comparison to the results of the 
single waveguide will be made. It is important to distin- 
guish between two cases: a) The monomer correlations, 
where we talk about the correlations between the fields 
within a single waveguide given by Eq. (^2j) and where 
coupling is still present, b) The limit of no coupling, 
where the spectra will behave as a single isolated waveg- 
uide. This limit is important since it allows us to com- 
pare with the results previously obtained by Collett and 



Walls ||, and henceforth this limit is referred to as the 
SHG monomer. Finally, we denote the spectral variances 
Vj^g k (u>) as the dimer correlations or variances. 

It was shown by Collett and Walls Q that in the 
SHG monomer without detuning very good squeezing 



in the fundamental quadrature —i(Aie 



A[e 



) is 



obtained when 7 is small, and conversely when 7 is 
large good squeezing in the second-harmonic quadrature 
— i(A 2 e _ * e2 — Ale 1 ® 2 ) is observed. These squeezing spec- 
tra were optimized by choosing a proper value of the 
quadrature phase and as it turns out 61 = 2 = 7r/2 
maximizes the squeezing in both cases (corresponding to 
the amplitude quadrature). For exactly this value the 
quadrature correlations coincide (to leading order) with 
the monomer intensity correlations (^lj) so the results 
of Ref . !| also predicts excellent noise suppression in the 
monomer photon number variances considered in this pa- 
per. Note that the choice 8 = n/2 only maximizes the 
squee zing when detuning is zero, as it was shown by Olsen 
et al. |42[| . 

Generally, the violation of the classical (or shot-noise) 
limit requires that the fluctuations diverge in a given 
observable of the fields. When this happens the spec- 
tral variance for this observable becomes large, and the 
canonical conjugate observable of the fields will in turn 
have a small variance as a consequence of the Heisen- 
berg uncertainty relation of canonical conjugate observ- 
ables. A typical situation where the fluctuations diverge 
is close to a transition from one stable state to another, 
and therefore violations of the shot-noise limit is nor- 
mally studied close to bifurcation points. In this paper 
we study the sum and the difference of the intensities of 
the fields, so a violation of the classical limit implies that 
sub-Poissonian statistics is observed and that the pho- 
tons at the photodetectors are anti-bunched; they arrive 
more regularly than if coherent beam intensities (which 
obey Poissonian statistics) were measured. The problem 
with the intensity observable is to find the conjugate ob- 
servable in which the fluctuations should become large 
when the intensity correlations violate the classical limit. 
Numerous attempts to create the most intuitive conju- 
gate observable, namely a phase operator, has not been 
entirely successful J3(J . On the other hand, in a photon- 
counting experiment it is exactly intensity correlations 
that are measured, making them a suitable choice for a 
direct experimental implementation. 
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The analytical and numerical results presented in the 
following display excellent mutual agreement. In order to 
achieve this it was necessary to have the time resolution 
of the two-time correlations low enough to describe the 
temporal variations, while simultaneously keeping upper 
limit of the correlation time (corresponding to the limits 
t — * ±00 of the analytical integral) long enough for the 
two-time correlation to become close to zero. Otherwise, 
the temporal Fourier transform of the correlations will 
give spectra that are in disagreement with the analytical 
results. Needless to say, this had to be checked for each 
case as the parameters were varied, but generally we used 
N = 512 or 1024 points with a resolution in the range 
At = 0.04 — 0.2 to calculate the two-time correlations 
C(t). Finally, the length of the simulations was around 
10 6 time units before the correlations converged to the 
degree shown in the following. 
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FIG. 4. Photon number spectra V^ B (uo) for the param- 
eters in Fig. I and Ji = 3.0, A = and E = 3.275, on 
the lower branch just before a bistable turning point. Lines 
show analytical results while points are numerical results. The 
shot-noise level is indicated with Csn. 



A. 7 small 

First we consider the case where 7 is small and good 
nonclassical correlations in the FH fields are expected 
(Va 1 (uj) — 7/[l + 7] close to self-pulsing transitions) ||. 
For 7 = 0.1 we observed the strongest violations of the 
quantum limit close to bistable turning points. An ex- 
ample is shown in Fig. || located in the bistable region of 
Fig. [| and the system is set on the lower branch of the 
bistable curve just before the right turning point. The 
dimer spectrum of the sum of the FH fields shows a near 
Lorentzian dip in the region of lu = that goes down 
to vj^ Bi (cu) ~ 0.2, implying strong twin-beam correla- 
tions, while the FH difference shows excess noise here. 
Taking 7 even smaller we were able to get V A ^ B (lu) very 
close to zero in the presence of bistable turning points, 
a behavior similar to the 7/(1 + 7) behavior observed in 
the SHG monomer close to self-pulsing transitions. The 



excellent correlations are only seen close to the bistable 
transition, taking e.g. E = 3.2 for the parameters in 

Fig. |], the minimum of the spectrum is Vj^ B (to) ~ 0.35. 
Returning to Fig. |, at w ~ 4 the FH sum spectrum 
again shows nonclassical correlations of approximately 
60% of the shot-noise limit. The frequency almost co- 
incides with the frequency of the Hopf instability that 
is suppressed here; the eigenvalue with the largest real 
part is the bistable one while the eigenvalues with more 
negative real parts have Im(A) ~ 3.6. 

Here it is relevant to mention that bistability is also 
present in the SHG monomer |3^,|^| (as Eq. ( |22] ) indi- 
cates this requires nonzero detunings with equal sign), 
and to the best of our knowledge nobody has here inves- 
tigated the quantum behavior. Let us write the detunings 
of the SHG monomer as A,-. Due to the invariance of the 



symmetric steady state solutions when Aj 



Ji 



we can obtain the same bistable state investigated in 
Fig. ^ in the SHG monomer by setting Ai = —3.0 and 
A2 = —1.0. The spectrum Va 1 (ui) displays here exactly 

the same behavior as Vj£ B (to) in Fig. ||, so also in the 
SHG monomer perfect anti-bunching behavior may be 
obtained in the small 7 limit. Generally, the dimer spec- 
tra vj$^ B . (lu) can be reproduced by the monomer spectra 



(to) when taking Aj 



Jj. This is not valid, 



however, close to a transition to asymmetric states, and 
also the spectra V a . b .(lo) have no equivalents in the no- 
coupling intensity correlations. 
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FIG. 5. Photon number spectra V^ B (co) for the param- 
eters in Fig. H Ji = 2.0, A = and E = 2.3. Lines show 
analytical results while points are numerical results. The 
shot-noise level is indicated with Csn- 

In the self-pulsing region of Fig. || it is possible to ob- 
tain good correlations if the system is set close to the 
bistable area. The spectra in Fig. are for a pump 
value where both the bistable and the self-pulsing eigen- 
values are of approximately the same strength, and the 
plot shows that the dimer spectrum of the FH differ- 
ence have strong noise suppression for nonzero w, that 
originates from the self-pulsing instability setting in at 







Esp = 4.7. Also the sum shows strong noise suppression 
now at w = 0, caused by the proximity of the bistable 
area which gives rise to an eigenvalue with Im(A) = that 
never has Re(A) > 0. The good correlations observed 
here are apparently a result of a competition between the 
bistable state and the emerging self-pulsing instability, 
that eventually dominates for higher pump levels. When 
the self-pulsing threshold is approached, the nonclassical 
behavior is less pronounced. This does not necessarily 
mean that nonclassical states are not present here, but 
probably that the intensity is the wrong observable in 
which to observe nonclassical correlations here. Such a 
case was observed in the SHG monomer, where it was 
found that detuning of the system makes the squeezing 
ellipse turn so the best squeezing is no longer observed 
in the amplitude quadrature E2fl, 
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FIG. 6. Photon number spectra for Ai = A2 = 1.1 and 



7 



0.1, Ji = 20.0, J 2 = 1 and E/E s: 



0.97. Lines 



show analytical results while points are numerical results. The 
shot-noise level is indicated with Csn- 

When detuning is introduced, the system can be- 
come asymmetrically unstable for certain couplings as 
was shown in Fig. |3j. Generally, the asymmetric in- 
stability shows sub-Poissonian twin-beam correlations in 
the dimer FH difference spectrum, which is especially 
pronounced when J\ 3> J2 where almost perfect anti- 
bunching was observed. In Fig. ^ the dimer spectra 
vj^g (uj) are shown for A = 1.1, J\ = 20 and J 2 = 1, 
taken close to the symmetric transition E sym , and inten- 
sity correlations until 8% of the shot-noise limit is seen in 
the FH difference correlations at nonzero u>. By carefully 
selecting the parameters we were even able to see corre- 
lations until 3% of the shot-noise level, which underlines 
that excellent nonclassical correlations are observed here. 
This result is quite robust; good sub-Poissonian correla- 
tions are observed also further below the transition as 
well as for considerably lower values of the FH coupling 
strength, while the peaked structure around lo = is 
quite sensitive to the pump level since it is not seen tak- 
ing the system even closer to E sym . We note from Fig. || 
the presence of the self-pulsing instability for the param- 
eters chosen, and even if quite good correlations are ob- 



served inside the large asymmetric area for J\ large, the 
best results are obtained close to the self-pulsing regions. 
Again this shows that two competing instabilities appear 
to give rise to strong nonclassical correlations. We note 
finally that the frequency u> ~ 11, where the best corre- 
lations in Fig. [| are observed, almost coincides with the 
frequency of the self-pulsing eigenvalue that is damped 
the most. Thus, paradoxically, in this case the least dom- 
inating eigenvalue is determining the frequency of the 
strong correlations. 
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FIG. 7. Photon number spectra on resonance and for 7 = 1 
and Ji = 4.0, J2 = 1.0 and E/Esp = 0.95. Lines show analyt- 
ical results while points are numerical results. The shot-noise 
level is indicated with Csn- 



B. 7 = 1 

Setting the loss rates to be identical, 7 = 1, the self- 
pulsing instability gives rise to strong nonclassical corre- 
lations all the way to the transition to the self-pulsing 
state, while the bistable transition displays only weak vi- 
olations. As an example of the self-pulsing correlations 
Fig. [7] displays selected spectra for the dimer correlations. 
The sum of the SH fields displays strong correlations at 
u = Q, which goes to 25% of the shot-noise limit when 
E — > -EsP: a result that can also be found in the SHG 
monomer by introducing equivalent dctunings. Also the 
dimer FH difference correlations arc strong; for nonzero uj 
correlations around 50% of the shot-noise limit are seen. 



C. 7 large 

For large 7 the SHG monomer predicts strong shot- 
noise violations in the SH spectrum at the self-pulsing 
threshold [bj (Va 2 (w) — (1 + 7) -1 close to self-pulsing 
transitions), and for the dimer similar levels of corre- 
lations can be obtained close to the self-pulsing transi- 
tion. As an example of the behavior Fig. H shows highly 
nonclassical twin-beam correlations in the dimer spec- 
trum for the sum of the SH intensities. For the selected 
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parameters we observe also sub-Poissonian behavior in 
Vj^g (u), which, in contrast to Vj£ B (cv), cannot be ob- 
served in the SHG monomer with equivalent detunings. 
Although asymmetric areas were found for nonzero de- 
tunings, only weak sub-Poissonian correlations were ob- 
served there in the difference of the SH fields (the sum 
correlations still show strong sub-Poissonian correlations 
here, as expected from the SHG monomer predictions). 
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FIG. 8. Photon number spectra on resonance and for 
7 = 10 and Jj = 6.0, J 2 = 2.0, and E/E SP = 0.95. Lines 
show analytical results while points are numerical results. The 
shot-noise level is indicated with Csn- 



VI. SUMMARY AND DISCUSSION 

In this paper we have proposed a model we denote the 
quantum optical dimer for studying the effects of a simple 
mode coupling in a cavity. The model consisted of two 
%( 2 ) nonlinear waveguides in a cavity, with coupling be- 
tween them from evanescent overlapping waves that was 
assumed weak and linear. We chose to restrict ourselves 
to investigating nonclassical correlations, and hence de- 
rived the nonlinear quantum equations of evolution for 
the system, resulting in a set of stochastic Langevin equa- 
tions. 

Using a linearized analysis we showed that the system 
for low pump levels allowed a symmetric state to be sta- 
ble, where both waveguides are in the same stable state. 
Depending on the system parameters this state destabi- 
lized into a self-pulsing state, where temporal oscillations 
are observed, or bistable solutions in the steady states oc- 
curred. For some parameters, the symmetric steady state 
lost stability in favor of an asymmetric steady state, in 
which the two waveguides has different steady state so- 
lutions. 

We investigated the effects of the quantum noise 
present in the system by calculating two-time intensity 
correlation spectra of the output fields. It was shown 
that sub-Poissonian correlations were present in the sys- 
tem, especially when considering the sum or difference of 
the field intensities from each waveguide implying strong 
twin-beam correlations. This nonclassical anti-bunching 



effect is a true manifestation of a quantum behavior and 
was observed mainly in three cases: 

• Close to bistable turning points the strongest vio- 
lations of the classical limit were observed in the 
spectrum at u = 0, corresponding to correlations 
at infinite time. In the limit of 7 <C 1 perfect noise 
suppression could be obtained in the sum of the FH 
intensities, resulting in perfect twin-beam behav- 
ior. This was also observed in the single waveguide 
model, when equivalent detunings were introduced. 

• Close to threshold for self-pulsing behavior strong 
twin-beam correlations were observed both at w = 
and also at values of to coinciding with the oscil- 
lation frequency of the emerging instability. Thus, 
the correlations here anticipate the behavior above 
the threshold analogously to the idea of a quantum 
image [Q, where the spatial modulations are en- 
coded in the correlations below threshold while the 
average intensity remains homogeneous. 

• Excellent sub-Poissonian correlations were ob- 
served close to transitions from symmetric to asym- 
metric steady states, which is a wholly unique tran- 
sition to the dimer. This instability is closely re- 
lated to the near-field of a modulationally unsta- 
ble system in presence of diffraction; the dimer 
sites could be thought of as neighboring near-field 
pixels. Variances down to 3% of the shot-noise 
level were observed in the the difference of the 
FH fields, implying nearly perfect twin-beam be- 
havior. The correlations were particularly strong 
when the FH coupling strength was much larger 
than the SH and when the FH loss rate was much 
larger than the SH loss rate. These results indicate 
that strong near-field correlations can be observed 
in SHG with diffraction close to a modulational in- 
stability, which to our knowledge has not been in- 
vestigated yet. 

It is worth noting that the twin-beam correlations re- 
ported here were all originating from the dimer coupling 
across the waveguides; while each field was created in- 
dividually from the nonlinear interaction in the corre- 
sponding waveguide, the coupling between the waveg- 
uides gave rise to the strong nonclassical twin-beam cor- 
relations. Hence, the nonclassical correlations arise not 
because the photons are twins, but rather because they 
are "brothers". 

Common for all these cases was that distinctively 
strong nonclassical correlations were observed in param- 
eter regimes where two types of instabilities were com- 
peting. This was observed in self-pulsing areas close to 
bistable and asymmetric regimes, and also in asymmet- 
ric areas close to bistable and self-pulsing regimes. These 
enhanced correlations from competing instabilities has to 
the best of our knowledge not been reported before. 
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We showed that the relative input mirror loss rate 
7 = 72 /71 between the SH and FH fields had a strong in- 
fluence on the sub-Poissonian behavior, as was previously 
shown by Collett and Walls M in the SHG monomer. For 
small 7 the strongest nonclassical states were mainly ob- 
served in the FH fields while for large 7 they were mainly 
observed in the SH fields. Since the photon life times in 
the cavity are inversely proportional to the loss rates of 
the input mirror, the time spent in the cavity is decisive 
for the level of nonclassical correlations of the output 
fields; the field with the shortest time spent in the cavity 
displays the strongest nonclassical correlations. The fact 
that a long interaction time tends to destroy the nonclas- 
sical correlations have also been observed in propagation 
setups, specifically Olsen et al. showed that in propa- 
gation SHG the presence of quantum fluctuations caused 
a dramatic revival of the FH after a certain propagation 
length, causing the variance to go above shot-noise level. 

We stress that the results presented here were very 
robust to changes in the parameters. Thus, large param- 
eters areas exist where strong nonclassical behavior can 
be seen. This is especially important to stress for the cou- 
pling parameters, since they are not so easily controlled 
experimentally. 

We investigated only the symmetric state of the sys- 
tem, and a future study of the stability, dynamics and 
nonclassical properties of the asymmetric states is rele- 
vant. In this context we should also mention that the 
coupling between the waveguides chosen here was con- 
servative (imaginary coupling in the Langevin equations) 
while it could also have been dissipative (real) or a com- 
bination (complex). This highly depends on the actual 
setup, and a future study should include the possibility 
of a general complex coupling. 
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APPENDIX A: THE QUANTUM-TO-CLASSICAL 
CORRESPONDENCE 

In this section we show how the master equation (|^) 
is converted into an equivalent partial differential equa- 
tion by expanding the density matrix in coherent states 
weighted by a quasi-probability distribution pq , p9] . In 
this distribution the operators are replaced by equivalent 
c- numbers, where the particular correspondence between 
these depends on the ordering of the operators. In the 



case where only up to second order derivatives appear in 
the corresponding partial differential equation, the equa- 
tion is on Fokker-Planck form allowing equivalent sets of 
stochastic Langevin equations to be found. 

We are now left with a choice of probability distribu- 
tion, be it either the P, Wigner or Q distribution giving 
normal, symmetric or anti-normal averages, respectively. 
The aim of this paper is to calculate two-time correla- 
tion spectra outside the cavity, and in order to do this 
most conveniently, the moments of the intracavity fields 
must be time and normally ordered, cf. the discussion 
in Sec. IV. Since the P- representation will immediately 
give the time and normally ordered averages needed, this 
is the favorable representation in this context. As will 
be explained later we use the Wigner representation for 
the numerical simulations, and therefore we now provide 
a general way of deriving the QPD equations. 

The approach we use is to introduce a characteristic 
function 



X (z) = Tr[D(z)p], 



(Al) 



so x is the trace over a displacement operator D{z) act- 
ing on the density matrix [i.e. the expectation value of 
D(z)\. The choice of ordering now amounts to choosing 
the ordering of D(z). In the symmetric ordering 



D„(z) = e z 



(A2) 



where z is a complex number describing the amplitude of 
a coherent field, and A is a boson operator. The normally 
ordered displacement operator is 

D n (z) = e zAi e-** A , (A3) 

and the anti-normally ordered displacement operator is 

D a (z) = e- z * A e zA \ (A4) 



The QPD is now given as a Fourier transform of the char- 
acteristic function 



W{a) 



d 2 zx(z)e 



z a — za 



(A5) 



where the integration measure d 2 z means integration over 
the entire complex plane. From this relation an equiva- 
lence between operators and c-numbers has been estab- 
lished as A «-> a and <-> a*. The c-number averages 
may now be calculated as, e.g., 



{a a) 



d 2 aW(a)a*a, 



and obviously since the operator o rder ing determines the 
specific form of W(a) from Eq. (A5), this in turn im- 
plies that the c-number averages also are influenced by 
the choice of ordering. 
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From differentiating Eq. (Al) with respect to time we 



get 



dx(z) 
dt 



Tr D(z) 



dp 

~dt 



(A6) 



with dtp being governed by the master equation (|^), and 
we have used that in the Schrodinger picture the opera- 
tors [here D{z)\ are independent of time. The approach 
is now to differentiate D(z) with respect to e.g. z and 
rearrange to get 



d_ 

dz 



D s (z) 



(A7) 



The righ t hand side of Eq. ( |A6| ) is evaluated us ing 
Eq. (A7) and the similar other expressions. Eq . ( |A6|) is 
then Fourier transformed according to Eq. (A5), assum- 
ing the characteristic function is well behaved, to finally 
give the equation governing the time evolution of W(a). 



Choosing the no rma lly ordered displacement opera- 
tor given by Eq. (A3) the equation for the Glauber- 
Sudarshan P-rcprcscntation is derived, which is on 
Fokker-Planck form. However, due to problems with 
negative diffusion in quantum optics the generalized P- 
distributions |^,^| are normally used instead, where 
the problems are surpassed by doubling the phase space. 
We will use the positive P-representation, which can 



be derived by replacing all a* — > a] and j3* — > /?j in 
the Fokker-Planck equation of the Glauber-Sudarshan 
P-representation. This means that aj is now an inde- 
pendent complex quantity instead of being the complex 
conjugate of ctj. The Fokker-Planck equation using the 
positive P-representation corresponding to the master 
equation (||) is then 



aw n (x) d 

dt | dai 



0*1(71 — iSi) + iJxfii — na\a 2 — £ p 
a l(7l + *^i) — iJiPi — Kai«2 ~~ £p, a 



012(72 - 182) + 1J2P2 + -jOL\ 
a\(p(2 + iS 2 ) - U2P2 + 7^a\ 



/?l(7l _ »<*l) + iJiUi - K/3{P 2 ~ £p,b 
/3j(7i + i8i) - iJia\ - K,f}\fr 2 ~ £p,b 



^2(72 - iS 2 ) + i J2U2 + 
4(72 + iS 2 ) - Uia\ - |/3| 2 



Jtt2 



2 "2 



+^A + -^a4]jw„(x), 
d Pi d(3\ J J 



(A8) 



where the c-number equivalents of the operators are 
{A V A]} {a v a]} and {P,,Pj} «-> and x 

represents the c-number states 

x = {ai, a\,a 2 ,a\, Pi,f3{, f3 2 ,p\}. 



The numerical simulation of the positive P- 
representation has been reported as very difficult, mainly 
due to divergent trajectories 48 1, cf. the discussion in 
App. [b|. Therefore we choose to use the Wigner rep- 
resentation for the numerical simulations, obtained by 
using the symmetric displacement operator (A2). The 
time evolution of the Wigner distribution is governed by 



dW 3 (x) j d 



dt 



dai 



[(71 - iSi)ai - na\a 2 + iJifii - £ P)i 



+ 77^- [(72 - 162)012 + ^-a\ + iJ 2 (3 2 ] 

oa 2 i 

d 

+ ~^7T [(7i - ^i)/5i - + iJ\ol\ - £ p : 

dpi 

+ flo"[(7a " iS 2)p2 + *0£ + U 2 a 2 
dp 2 



2 

d 2 



\da\da\ 




( d 2 


d 2 


\da 2 da2 


dfodfc 


' d 3 


d 3 > 



4 \da 2 da* 2 d(3 2 d(3* 2 



c.c. ^(x), (A9) 
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where the c-number equivalents of the operators are 
{A*, i}} «-> {a 3 ,a*} and {B^b]} {,%/?*} and 

x = {ai, 0%, a 2 , oil, 0i, 01, /3 2 , (3%}. 



Due to the third order derivatives Eq. (|A9| ) is not on 
Fokker-Planck form, a problem we address in App. 
Note that the +c.c. term (denoting the complex conju- 
gate) at the end applies to the entire equation. 

The connection from the QPD equations to the 
stochastic Langevin equations can be made if the QPD 
equation is on Fokker-Planck form, which for a system 
with m c-number states Xj is p9| 



cW(x) 

Ft 



2 dxjdxu 

— 



D 



(x)}w(x), (AlO) 



Using Ito rules for stochastic integration the equivalent 
set of Langevin equations is 



dx ' 

-^=A j (x)+w j (t), 



(All) 



where utj(t) are Gaussian white noise terms, delta corre- 
lated in time according to the diffusion matrix D 



(w j {t)w k (t , ))=D 3k {x)5(t-t r ). 



(A12) 



We note that if D depends on x the noise is labeled mul- 
tiplicative which is the case for the positive-P Eq. ( AS ), 
otherwise it is additive as it is for the Wigner Eq. (A9). 



APPENDIX B: NUMERICAL SIMULATIONS 

The choice of using the Wigner representation for the 
numerical simulations is not immediately apparent, since 
it involves an approximation that is not necessary if the 
P or the Q-representation are used. The advantage of 
the truncated Wigner Langevin equations (0) is that the 
noise is additive, as opposed to the the multiplicative 
noise of the Q" re P resei rtation (where the noise for the 
quantum SHG model poses serious limits on the param- 
eter space and the P-representation. For the P- 
representation we are forced to use the generalized rep- 
resentations in order to avoid negative diffusion in the 
Fokker-Planck equation. Since this choice implies dou- 
bling of the phase space, the c-numbers ctj and ctj are no 
longer each others complex conjugate (only on average) 
and the respective noise terms are not correlated to each 
other, cf. Eq. ([To]). This may lead to divergent trajecto- 
ries where the convergence is extremely slow, and is the 
major reason for us avoiding a numerical implementation 
of the positive-P equations. The Wigner equations, on 
the other hand, have no problems in this direction. 



The drawbacks to using the Wigner equations are first 
of all that we have to ne glect the third order terms of the 
Wigner QPD equation (A9) to get it on Fokke r-Planck 
form so the equivalent Langevin equations ( All ) may be 
obtained. It is uncertain what the implications of this 
approximation are, however in many cases no major dif- 
ferences have been observed between simulations of the 
truncated Wigner equations compared to exact positive- 
P or Q equations (4^jT^|. On the other hand in Ref. 
fl32t so-called quantum jump processes in the degener- 
ate OPO above threshold are shown to produce signifi- 
cant differences between the truncated Wigner and the 
positive P-representation. In our case the third order 
terms are 0(k 4 ) while the other terms are 0(k 2 ) or lower. 
And because of the weak nonlinear coupling the effect of 
the third order terms is weak, justifying the truncation. 
Another drawback to the Wigner representation is that 
the intracavity averages are symmetrically ordered, and 
these cannot be rewritten to time and normal ordering 
since the intracavity commutator relations are not known 
for t t' (only the output fields have well defined cor- 
relations here). This means that in order to compute 
the output fields at a given time t we are forced to use 
Eq. (|2^) , and here the Gaussian white noise part of the 
input field is an ill-defined instantaneous quantity. The 
output fields of the numerics are calculated by using the 
fact that although the derivative and instantaneous value 
of a stochastic term arc ill-defined, the integral is well- 
defined. By integrating over a time window (which we 
denote At) and calculating the average, as described in 
Ref. |49), we may obtain the output fields from Eq. (|2^). 

We use the Heun method |k| to numerically solve the 
Langevin equations for the intracavity fields and to eval- 
uate the output fields, and a random number generator 
pl| for generating the Gaussian noise terms. The time 
step was set to At = 0.001 and checked to be stable. The 
size of the time window for calculating the output fields 
was varied between At = 40 At — 200Ai according to the 
resolution needed for the individual spectra. Finally, we 
set the noise strength parameter n s = 10 8 , which is a 
typical value for the cavity configuration considered here 

The averages calculated using the output fields of the 
numerics correspond to symmetrically ordered averages 
since they are calculated from the Wigner Langevin equa- 
tions. In order to relate these averages to the normally 
ordered averages of the spectra in Sec. |y|, the output 
commutator relations (HJ) are used to rewrite the out- 
put correlations. The classical steady states of the output 
fields are found from the average of Eq. ( p3|) as 




2 7l ^ 1 +£ p / v /2^r 
2^F 2 , T = A,B, 



(Bla) 
(Bib) 



by taking the input fluctuation to be zero on average. 
Assuming that the output fields are fluctuating around 
the output steady states 



j, out i 



(B2) 
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we introduce the photon number fluctuation operator for 
waveguide A as 



<^out> 



A 



-j.outAA] out (t) 



A*°ut A ^\°utW, 



(B3) 



with subscript s to indicate that the average is symmetric 
and we have neglected higher order terms in the fluctu- 
ations. Using this expression, the two-time correlation 
function (p7|) is with a symmetric ordering of the opera- 
tors to leading order 



C 



£> (r) ~ <AJV/ out (i) ± AiStf 



j .out V 7 — — 1 ' ' /, ,ut (0' 

.4 /. , _\ i A XV-B 



AA^ out (t + r) ± AiV£ out (* + r)) s . (B4) 



The symmetric c-number correlations of the output 
field fluctuations are now calculated in the numerical 
simulations of the dimensionless Wigner Langevin equa- 
tions (0) as (Aw s (0), [Aw s (t)] t ) s where 

Aw s = [AA 1 ,AA* 1 ,AA 2 ,AAl,AB 1 ,ABl,AB 2 , AB*] T . 

The spectral matrix S s (w) of fluctuations is now straight- 
forwardly given by the Fou rier transform of these correla- 
tions, and the correlations (B4) may now be calculated in 
the same manner as shown in Sec. IV with the shot-noise 
level 



C SN = (\A 



out ,j I 



IS, 



out,/c| 



2 K 



using that the normalization of the output fields ar e th e 
same as the one taken for the input fields in Eq. (11c). 
Note that the shot noise level, here expressed in the sym- 
metric averages in the Wigner Langevin equation, is iden- 
tical to the shot-noise level expressed in averages from the 
positive-P equations (35), since |-4 ut,j| 2 = ^Ijlj- This 



is due to the approximation made in Eq. (B3). 

In the analytical treatment in Sec. IV we used that in 
the spectral matrix S(w) certain symmetries are present, 
so in fact only approximately one third of the 64 correla- 
tions were needed to obtain the results presented there. 
In a numerical simulation this is only approximately valid 
in the limits of long integration times and large correla- 
tion times. Much better results are obtained faster if the 
spectra are calculated directly from the full 8x8 matrix 
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